\section{Divergence Free Field In Fourier Domain}
A brief proof for divergence free vector field's property in Fourier domain: Its direction in the Fourier domain is always perpendicular with the wavenumber direction.

Suppose we have a 2D vector field $\textbf{u} = [u_x, u_y]^T$ and the Inverse Fourier Transform of each component is
\begin{equation}
\begin{split}
u_x(x, y) &= \frac{1}{N\times M}\sum_{p=0}^{N}\sum_{q=0}^{M}e^{i(w_px+w_qy)}\hat{u}_x(p,q)\\
u_y(x, y) &= \frac{1}{N\times M}\sum_{p=0}^{N}\sum_{q=0}^{M}e^{i(w_px+w_qy)}\hat{u}_y(p,q)\\
\end{split}
\end{equation}

Here N and M are the resolution of the discrete Fourier domain. $\hat{\textbf{u}}=[\hat{u}_x,\hat{u}_y]^T$ is \textbf{u} in Fourier domain, which is also a 2D vector field. In order to get a divergence free field, we want \textbf{u} to satisfy $\nabla\cdot\textbf{u}=0$ or $\frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} = 0$. In this sense, we have the equivalent form in Fourier domain
\begin{equation}
\begin{split}
\frac{\partial u_x}{\partial x} &= \frac{1}{N\times M}\sum_{p=0}^{N}\sum_{q=0}^{M}iw_pe^{i(w_px+w_qy)}\hat{u}_x(p,q)\\
\frac{\partial u_y}{\partial y} &= \frac{1}{N\times M}\sum_{p=0}^{N}\sum_{q=0}^{M}iw_qe^{i(w_px+w_qy)}\hat{u}_y(p,q)\\
iw_pe^{i(w_px+w_qy)}\hat{u}_x(p,q) &= w_p\hat{u}_x(p,q)e^{i(w_px+w_qy+\frac{\pi}{2})}\\
iw_qe^{i(w_px+w_qy)}\hat{u}_y(p,q) &= w_q\hat{u}_y(p,q)e^{i(w_px+w_qy+\frac{\pi}{2})}
\end{split}
\end{equation}

To satisfy divergence free condition, we should have

\begin{equation}
w_p\hat{u}_x(p,q) + w_q\hat{u}_y(p,q) = 0
\end{equation}

because the Fourier basis functions are orthonogal basis. Eqn. 3 is simply the dot product between $[w_p, w_q]^T$ and $\hat{\textbf{u}}$, which implies they are perpendicular with each other. And $[w_p, w_q]^T$ is $[p, q]^T$ scaled by some scalar. Thus we get to the conclusion.